#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use chromosomes;
use linewrap;

use db_parameters;
use ortho_species;

print <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    paralogues functional summary

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

#
#	Connect to genome database
#
my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);
$dbh->{RaiseError} = 0;
#my @db_param = get_db_parameters();



#_________________________________________________________________________________________

#	get_clusters_above_size


#	Retrieve clusters in order of cluster sizes

#_________________________________________________________________________________________
sub get_clusters_above_size($$$)
{
	my ($dbh, $cluster_size, $protocol_id) = @_;
	my @clusters;
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
        ortholog_id  
    FROM 
        orthologs.ortholog_sets 
    WHERE
        ortholog_type ~ 'to' AND
        protocol_id = $protocol_id
    GROUP BY ortholog_id
	HAVING count(*) > $cluster_size;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for (@$data)
	{
		push(@clusters, $_->[0]);
	}
	return @clusters;
}


#_________________________________________________________________________________________
#
#	retrieve_interpro_data

#		Retrieve interpro names and descriptions for all matches

#_________________________________________________________________________________________
sub retrieve_interpro_data($\%)
{
	my ($dbh, $interpro_data) = @_;
	print STDERR "\tRetrieve interpro names and descriptions for all matches...\n";
    my $sql_cmd = <<"PL/SQLCMD";
    SELECT
        interpro_ac,
        interpro_name,
        interpro_description
    FROM
        domains.interpro_data;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
    {
        $interpro_data->{$line->[0]} = [$line->[1], $line->[2]];
    }
}


#_________________________________________________________________________________________
#
#	retrieve_taxon_prot_interpro_matches

#_________________________________________________________________________________________
sub retrieve_taxon_prot_interpro_matches($\%)
{
	my ($dbh, $interpro_per_protocol_per_ortholog) = @_;
	print STDERR "\tObtain interpro matches for each ortholog set:\n";
	print STDERR "\tSort interpro matches from $ortho_name1, $ortho_name2...\n";
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
        DISTINCT ON (protocol_id, ortholog_id, interpro_ac) 
        protocol_id, ortholog_id, interpro_ac
    FROM
        taxon.prot_interpro_matches NATURAL JOIN
        orthologs.ortholog_sets a
    WHERE
        species = '$ortho_name1' OR
        species = '$ortho_name2';
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		push(@{$interpro_per_protocol_per_ortholog->{$line->[0]}{$line->[1]}}, $line->[2]);
	}
}


#_________________________________________________________________________________________
#
#	retrieve_go_matches

#_________________________________________________________________________________________
sub retrieve_go_matches($\%)
{
	my ($dbh, $go_per_protocol_per_ortholog) = @_;
	print STDERR "\tObtain GO term matches for each ortholog set:\n";
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
            DISTINCT ON (protocol_id, ortholog_id, GO_id)
            protocol_id,
            ortholog_id,
            GO_id
        FROM 
            taxon.gene_go_matches gtpg
        NATURAL JOIN 
            (
                SELECT
                    GO_id
                FROM
                    orthologs.significant_go_terms
                WHERE 
                    category = 'duplicate' AND 
                    over_representation = TRUE
            ) AS a 
        JOIN
            orthologs.ortholog_sets USING (gene_id);
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		push(@{$go_per_protocol_per_ortholog->{$line->[0]}{$line->[1]}}, $line->[2]);
	}
}


#_________________________________________________________________________________________
#
#	retrieve_go_matches_per_geneid

#_________________________________________________________________________________________
sub retrieve_go_matches_per_gene_id($\%)
{
	my ($dbh, $go_per_gene_id) = @_;
	print STDERR "\tObtain GO term matches for each gene_id:\n";
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
            gene_id,
            GO_id
        FROM 
            taxon.gene_go_matches gtpg NATURAL JOIN 
            (
                SELECT
                    GO_id
                FROM
                    orthologs.significant_go_terms
                WHERE 
                    category = 'duplicate' AND 
                    over_representation = TRUE
            ) AS a;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		push(@{$go_per_gene_id->{$line->[0]}}, $line->[1]);
	}
}


#_________________________________________________________________________________________
#
#	retrieve_go_description

#_________________________________________________________________________________________
sub retrieve_go_description($\%)
{
	my ($dbh, $go_descriptions) = @_;
	print STDERR "\tRetrieve GO descriptions for all matches...\n";
    my $sql_cmd = <<"PL/SQLCMD";
    SELECT
        distinct on (go_id)
        go_id,
        description
    FROM
        orthologs.significant_go_terms NATURAL JOIN 
        go_term_descriptions gtd;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
    {
        $go_descriptions->{$line->[0]} = $line->[1];
    }
}






#_________________________________________________________________________________________

#	print_go_data

#_________________________________________________________________________________________
sub print_go_data($$$\%\%\%$)
{
	my ($dbh, $ortholog_id, $printH, $go_descriptions, $go_per_ortholog, $overrepresented_GOids, $species) = @_;
	return unless $go_per_ortholog->{$ortholog_id};
	my @matches = @{$go_per_ortholog->{$ortholog_id}};
	return unless @matches;

	print $printH "    GO matches for this cluster:\n";
	for my $go_id(@matches)
	{
		my $over = "";
		if (exists $overrepresented_GOids->{$go_id}{$species})
		{
			$over = ($overrepresented_GOids->{$go_id}{$species}[1]) ? "++over++" : "__over__";
		}
		die "[$go_id]\n" unless exists $go_descriptions->{$go_id};
		die "[$go_id]\n" unless defined $go_id;
		die "[$go_id]\n" unless defined $go_descriptions->{$go_id};
		printf  $printH "      %-15s %8s %-20s\n", "GO:$go_id", $over, $go_descriptions->{$go_id};
	}
	print  $printH  "\n";
}


#_________________________________________________________________________________________
#
#
#	retrieve_overrepresented_GO_ids
#
#_________________________________________________________________________________________
sub retrieve_overrepresented_GO_ids($\%)
{
	print STDERR "\tGet significantly overrepresented GO IDs...\n";
	my ($dbh, $overrepresented_GO) = @_;
	my $sql_cmd = <<'PL/SQLCMD';
    SELECT
        protocol_id,
		GO_id,
        species,
        CASE WHEN go_slim THEN 1 ELSE 0 END 
	FROM 
		orthologs.significant_go_terms
    WHERE
        category = 'duplicate' AND 
        significance < 0.0001 AND 
        over_representation = TRUE
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";

	# save GO identifiers per GO category per gene id
	for my $line (@$data)
	{
		++$overrepresented_GO->{$line->[0]}{$line->[1]}{$line->[2]}[$line->[3]];
	}

}



		

#_________________________________________________________________________________________

#	get_interpro_data

#_________________________________________________________________________________________
sub get_interpro_data($$\%\%\$)
{
	my ($dbh, $ortholog_id, $interpro_data, $interpro_per_ortholog, $interpro_hits) = @_;
	return unless $interpro_per_ortholog->{$ortholog_id};

	my @matching_descriptions;
	my @matches = @{$interpro_per_ortholog->{$ortholog_id}};
	for my $interpro_ac(@matches)
	{
		push(@matching_descriptions,  @{$interpro_data->{$interpro_ac}});
	}
	$$interpro_hits = join ("/", (@matching_descriptions));
}





#_________________________________________________________________________________________

#	gene loci per cluster

#		clusters are grouped by ortholog_id, species and adjacent clusters

#_________________________________________________________________________________________
sub get_per_cluster_loci($\%$$)
{
	my ($dbh, $cluster_loci, $species, $protocol_id) = @_;
	print STDERR ("\tRetrieving clusters for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		  ortholog_id,
		  local_cluster_id,
		  gene_id,
		  prot_id,
		  gene_position,
		  chromosome
	FROM
        orthologs.ortholog_sets NATURAL JOIN
		orthologs.local_paralogs NATURAL JOIN
		taxon.gene_loci NATURAL JOIN 
        orthologs.real_genes
    WHERE
        protocol_id = $protocol_id AND 
        ortholog_type ~ 'to' AND
        species = '$species'
	ORDER BY ortholog_id, local_cluster_id, gene_position;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		# {ortholog_id}{species}{local_cluster_id} = [gene_id, prot_id, gene_position, chromosome]
		push(@{$cluster_loci->{$line->[0]}{$species}{$line->[1]}}, [$line->[2],$line->[3],$line->[4],$line->[5]]);
	}
}

#_________________________________________________________________________________________

#	get_ortholog_cluster_ds

#_________________________________________________________________________________________
sub get_ortholog_cluster_ds($\%$)
{
	my ($dbh, $cluster_ds, $protocol_id) = @_;
	print STDERR ("\tRetrieving dS for each cluster...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
        ortholog_id,
        dS
    FROM
    	orthologs.ortholog_sets
    WHERE
        ortholog_type ~ 'to' AND
        protocol_id = $protocol_id
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	my %ortholog_id_to_ds;
	for my $line (@$data)
	{
		push(@{$ortholog_id_to_ds{$line->[0]}}, $line->[1]);
	}
	for my $ortholog_id(keys %ortholog_id_to_ds)
	{
		my @ds = sort {$a <=> $b} @{$ortholog_id_to_ds{$ortholog_id}};
		my $mid = int ((scalar @ds + 1) / 2);
		$cluster_ds->{$ortholog_id} = $ds[$mid];
	}
}



#_________________________________________________________________________________________

#	gene description

#_________________________________________________________________________________________
sub get_gene_descriptions($\%$)
{
	my ($dbh, $gene_descriptions, $species) = @_;
	print STDERR ("\tRetrieving gene descriptions for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		  gene_id,
		  description
    	FROM
    		taxon.gene_descriptions NATURAL join
            ens_id
        WHERE 
            species = '$species';    
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$gene_descriptions->{$line->[0]} =$line->[1];
	}
}

#_________________________________________________________________________________________

#	protein exon count

#_________________________________________________________________________________________
sub get_prot_exon_counts($\%$)
{
	my ($dbh, $exon_counts, $species) = @_;
	print STDERR ("\tRetrieving prot exon_counts for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
    		  prot_id,
    		  real_exon_count,
              disruptions_count
    	FROM
    		taxon.prot_exon_stats NATURAL join
            ens_id
        WHERE 
            species = '$species';
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$exon_counts->{$line->[0]} =[$line->[1], ($line->[2] ? $line->[2] : '')];
	}
}


#_________________________________________________________________________________________

#	protein exon count

#_________________________________________________________________________________________
sub get_prot_seqlen($\%)
{
	my ($dbh, $prot_len) = @_;
	print STDERR ("\tRetrieving prot lengths (residues) ...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		  prot_id,
		  len
	FROM
		orthologs.ortholog_sets NATURAL JOIN 
        accessions NATURAL join
        sequences;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$prot_len->{$line->[0]} = $line->[1];
	}
}


#_________________________________________________________________________________________

#	gene exon counts

#_________________________________________________________________________________________
sub get_gene_exon_counts($\%$)
{
	my ($dbh, $exon_counts, $species) = @_;
	print STDERR ("\tRetrieving gene exon_counts for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
            gene_id,
            exon_count
        FROM
            taxon.gene_coding_exon_counts NATURAL JOIN 
            ens_id
        WHERE 
            species = '$species';
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$exon_counts->{$line->[0]} =$line->[1];
	}
}




#_________________________________________________________________________________________
#
#	get_cluster_sizes

#	Retrieve clusters in order of cluster sizes

#_________________________________________________________________________________________
sub get_cluster_sizes($\%\@$$$)
{
	my ($dbh, $cluster_sizes, $clusters, $species1, $species2, $protocol_id) = @_;
	print STDERR ("\tRetrieving the list of clusters in descending order of $species1 cluster size...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
        ortholog_id, 
        sp1.paralog_count, 
        sp2.paralog_count
    FROM
        orthologs.paralog_stats sp1 join
        orthologs.paralog_stats sp2  USING (ortholog_id, protocol_id)
    WHERE
		sp1.paralog_count >= 1 AND
        sp1.species = '$species1' and
        sp2.species = '$species2' and
        protocol_id = $protocol_id
    ORDER BY 
        sp1.paralog_count desc;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$cluster_sizes->{$line->[0]} = [$line->[1], $line->[2]];
		push(@$clusters, $line->[0]);
	}
}

#_________________________________________________________________________________________
#
#	print cluster sizes

#_________________________________________________________________________________________
sub print_cluster_sizes(*\@\%\%\@)
{
	my ($printH, $ortholog_ids, $cluster_sizes, $cluster_ds, $species) = @_;
	print $printH "Cluster sizes for clusters > 6:\n";
	printf $printH "      %10s %11s %5s\n", "$species->[0] size", "$species->[1] size", "dS";
	printf $printH "      %10s %11s %5s\n", "----------", "---------", "-----";
	for my $ortholog_id(@$ortholog_ids)
	{
		printf $printH "      %10s %11s %5s\n",
					$cluster_sizes->{$ortholog_id}->[0],
					$cluster_sizes->{$ortholog_id}->[1],
					$cluster_ds->{$ortholog_id}
			if ($cluster_sizes->{$ortholog_id}->[0] + $cluster_sizes->{$ortholog_id}->[1]  > 6);
	}
	print $printH "\n\n\n";

}



#_________________________________________________________________________________________

#	print all summary info

#_________________________________________________________________________________________
sub print_summary(*\@\%\%\%\%\%\%\%\%\%\%\%\%\@)
{
	my ($printH,  $ortholog_ids, $cluster_sizes, $exon_counts, $gene_descriptions, 
		$cluster_loci, $cluster_ds, 
		$go_descriptions, $go_per_ortholog,
		$overrepresented_GOids,
		$interpro_data, $interpro_per_ortholog, $go_per_gene_id, $prot_len, $species) = @_;
	my $file_num = 0;
	print $printH join ("\t",  
						"ortholog_id",
						$species->[0],
						$species->[1],
						"description"), "\n";
	for my $ortholog_id(@$ortholog_ids)
	{
		last if ($cluster_sizes->{$ortholog_id}->[0] < 3);
		my $interpro_str = "";
		get_interpro_data($dbh, $ortholog_id, %$interpro_data, 
						  %$interpro_per_ortholog, $interpro_str);

		my %gene_descriptions;
		my %loci;
		for my $species($ortho_name1, $ortho_name2)
		{
			for my $sub_cluster_id(keys %{$cluster_loci->{$ortholog_id}{$species}})
			{
				for (@{$cluster_loci->{$ortholog_id}{$species}{$sub_cluster_id}})
				{
					# {ortholog_id}{species}{local_cluster_id} = [gene_id, prot_id, gene_position, chromosome]
					my ($gene_id, undef, undef, $chrm) =  @$_;
					++$loci{$species}{$chrm};
					next unless exists $gene_descriptions->{$gene_id};
					my $desc = $gene_descriptions->{$gene_id};
					$desc =~ s/\[.*\]//;
					$desc = lc $desc;
					++$gene_descriptions{$desc};
				}
			}
		}
		my $locus1 = join (",",  keys %{$loci{$ortho_name1}});
		my $locus2 = join (",",  keys %{$loci{$ortho_name2}});
		print $printH join "\t", ($ortholog_id, 
								  $cluster_ds->{$ortholog_id},
								  $locus1,
								  $locus2,
								  $cluster_sizes->{$ortholog_id}->[0],
								  $cluster_sizes->{$ortholog_id}->[1],
								  (join "/", keys %gene_descriptions), $interpro_str), "\n";
	}
	print STDERR "\tFinished summary information...\n";
}













#_________________________________________________________________________________________

#	get data

#_________________________________________________________________________________________

system ("mkdir -p $dir_orthologs_output/paralog_function");


# exon counts
my %exon_counts;
get_prot_exon_counts($dbh, %exon_counts, $ortho_name2);
get_prot_exon_counts($dbh, %exon_counts, $ortho_name1);
get_gene_exon_counts($dbh, %exon_counts, $ortho_name2);
get_gene_exon_counts($dbh, %exon_counts, $ortho_name1);

# gene descriptions
my %gene_descriptions;
get_gene_descriptions($dbh, %gene_descriptions, $ortho_name2);
get_gene_descriptions($dbh, %gene_descriptions, $ortho_name1);

# interpro data
my (%interpro_data, %interpro_per_protocol_per_ortholog);
retrieve_interpro_data($dbh, %interpro_data);
retrieve_taxon_prot_interpro_matches($dbh, %interpro_per_protocol_per_ortholog);


# go data
my %go_descriptions;
retrieve_go_description($dbh, %go_descriptions);
my %go_per_protocol_per_ortholog;
retrieve_go_matches($dbh, %go_per_protocol_per_ortholog);
my %go_per_gene_id;
retrieve_go_matches_per_gene_id($dbh, %go_per_gene_id);


#	Overrepresented GO identifiers for each protocol
#
my %overrepresented_GOids;
retrieve_overrepresented_GO_ids($dbh, %overrepresented_GOids);


# protein length data
my %prot_len;
get_prot_seqlen($dbh, %prot_len);

	
for my $protocol (get_orthologs_protocols())
{
	my ($protocol_name,  $protocol_id) = @$protocol;
	next unless $protocol_id == $curr_protocol_id;
	print STDERR "    Orthology protocol [$protocol_name]...\n";

	my %cluster_loci;
	get_per_cluster_loci($dbh, %cluster_loci, $ortho_name2, $protocol_id);
	get_per_cluster_loci($dbh, %cluster_loci, $ortho_name1, $protocol_id);
	
	
	
	my %cluster_ds;
	get_ortholog_cluster_ds($dbh, %cluster_ds, $protocol_id);
	
	#
	#	get cluster sizes
	#
	my %cluster_sizes1;
	my @cluster_ids1;
	get_cluster_sizes($dbh, %cluster_sizes1, @cluster_ids1, $ortho_name1, $ortho_name2, $protocol_id);
	
	my %cluster_sizes2;
	my @cluster_ids2;
	get_cluster_sizes($dbh, %cluster_sizes2, @cluster_ids2, $ortho_name2, $ortho_name1, $protocol_id);
	
	open SPECIES1_RESULTS,	">$dir_orthologs_output/paralog_function/orthologs.$ortho_name1.functional_summary" or die;
	open SPECIES2_RESULTS, 	">$dir_orthologs_output/paralog_function/orthologs.$ortho_name2.functional_summary" or die;
	my @species1 = ($ortho_name1, $ortho_name2);
	my @species2 = ($ortho_name2, $ortho_name1);
	#print_cluster_sizes(*SPECIES1_RESULTS, @cluster_ids1, %cluster_sizes1, %cluster_ds, @species1);
	#print_cluster_sizes(*SPECIES2_RESULTS,  @cluster_ids2, %cluster_sizes2, %cluster_ds, @species2);
	
	print STDERR "\tPrinting summary information by order of $ortho_name1...\n";
	print_summary(*SPECIES1_RESULTS, @cluster_ids1, %cluster_sizes1, %exon_counts, 
					%gene_descriptions, %cluster_loci, %cluster_ds, 
					%go_descriptions, %{$go_per_protocol_per_ortholog{$protocol_id}},
					%{$overrepresented_GOids{$protocol_id}},
					%interpro_data, %{$interpro_per_protocol_per_ortholog{$protocol_id}}, 
					%go_per_gene_id, %prot_len, @species1);
	print STDERR "\tPrinting summary information by order of $ortho_name2...\n";
	print_summary(*SPECIES2_RESULTS,  @cluster_ids2 , %cluster_sizes2 , %exon_counts, 
					%gene_descriptions, %cluster_loci, %cluster_ds, 
					%go_descriptions, %{$go_per_protocol_per_ortholog{$protocol_id}},
					%{$overrepresented_GOids{$protocol_id}},
					%interpro_data, %{$interpro_per_protocol_per_ortholog{$protocol_id}}, 
					%go_per_gene_id, %prot_len, @species2);
}


log_pipeline_stage($dbh, $curr_protocol_id, 46, 'table analyse paralogues', $start_time);
$dbh->disconnect();

